;Function to calculate surface area, given a particular grid
function gauge_flux_write_netcdf_areaGrid, lon, lat

  re=6367500.0  ;radius of Earth in m
  
  ;Assumes regular grid!
  dLon=abs(mean(Lon[1:n_elements(Lon)-1] - Lon[0:n_elements(Lon)-2]))*!pi/180.d
  dLat=abs(mean(Lat[1:n_elements(Lat)-1] - Lat[0:n_elements(Lat)-2]))*!pi/180.d
  Theta=double(!pi*(90.-Lat)/180.)
  
  AreaGrid=dblarr(n_elements(Lon), n_elements(Lat))
  
  For LatI=0L, n_elements(Lat)-1 do begin
    For LonI=0L, n_elements(Lon)-1 do begin
      if Theta[LatI] eq 0. or theta[LatI] eq !pi then begin
        Lat1=Theta[LatI] - dLat/2.
        Lat2=Theta[LatI] + dLat/2.
        AreaGrid[LonI, LatI]=(re^2.)*abs(cos(dlat/2.)-cos(0.))*dLon
      endif else begin
        Lat1=Theta[LatI] - dLat/2.
        Lat2=Theta[LatI] + dLat/2.
        AreaGrid[LonI, LatI]=((re^2.)*(cos(lat1)-cos(lat2))*dLon)
      endelse
    endfor
  endfor
  
  return, areaGrid
  
end

; docformat = 'rst'
;
;+
;
; :Purpose:
;   Create emissions netCDF file. 
;   This routine also checks that the grid is orientated correctly, 
;     and calculates and prints the global flux as a check that the units make sense.
;
; :Inputs:
;   lonInput: longitude
;   latInput: latitude
;   fluxInput: flux in kg/m2/s, dimensions [lon, lat, time]
;   time: array of Julian dates corresponding to each element in the time dimension of 'flux'
;   fileName: Name of netCDF file to be output
;   title: string detailing the contents of the file, used for the global "Title" attribute
;
; :Keywords:
;   author: author
;
; :Outputs:
;   Outputs 1 if successful
;
; :Example::
;   
;   IDL> lon=findgen(360) - 180.
;   IDL> lat=findgen(180) - 90.
;   IDL> time=julday(15, indgen(12)+1, 2000, 0)
;   IDL> flux=findgen(360, 180, 12)
;   IDL> success=gauge_flux_write_netcdf(lon, lat, time, flux, $
;     "/my/directory/N2O_flux_natsoil_2000_2001.nc", $
;     "N2O natural soil emissions from {citation} in kg/m2/s", $
;     author="Matt Rigby")
;   
; :History:
; 	Written by: Matt Rigby, University of Bristol, Jun 3, 2013
;
;-
function gauge_flux_write_netcdf, lonInput, latInput, time, fluxInput, fileName, title, $
    author=author

  compile_opt idl2
  on_error, 2

  lon=lonInput
  lat=latInput
  flux=fluxInput

  ;Check array dimensions
  if n_elements(flux[*, 0, 0]) ne n_elements(lon) then begin
    print, "NETCDF WRITE ERROR: Check flux longitude dimension (should be 1st dimension)"
    return, 0
  endif
  if n_elements(flux[0, *, 0]) ne n_elements(lat) then begin
    print, "NETCDF WRITE ERROR: Check flux latitude dimension (should be 2nd dimension)"
    return, 0
  endif
  if n_elements(flux[0, 0, *]) ne n_elements(time) then begin
    print, "NETCDF WRITE ERROR: Check flux time dimension (should be 3rd dimension)"
    return, 0
  endif

  ;Check that longitudes are ordered correctly
  whW=where(lon gt 180., countW, complement=whE)
  if countW gt 0 then begin
    print, "Longitudes should be ordered from -180 to 180, re-orientating"
    fluxSize=size(flux)
    fluxNew=fltarr(fluxSize[1:fluxSize[0]])
    fluxNew[0:countW-1, *, *]=flux[whW, *, *]
    fluxNew[countW:-1, *, *]=flux[whE, *, *]
    flux=fluxNew
    lon=[lon[whW] - 360., lon[whE]]
  endif

  ;Check that latitudes are ordered correctly
  if lat[1] lt lat[0] then begin
    print, "Latitudes should be ordered from -90 to 90, re-orientating"
    lat=reverse(lat)
    flux=reverse(flux, 2)
  endif

  ;Check global total flux
  area=gauge_flux_write_netcdf_areagrid(lon, lat)
  totalFlux=fltarr(n_elements(time))
  for ti=0, n_elements(time)-1 do begin
    totalFlux[ti]=total(flux[*, *, ti]*area)*365.25*24.*3600./1.e9
  endfor
  print, "Mean global flux: " + string(mean(totalFlux, /nan), format='(f8.4)') + $
    " Tg/yr. If this looks wrong, check units"

  ;Convert time
  timeType=size(time, /type)
  case 1 of    
    ;Julian dates
    (timeType eq 4 or timeType eq 5): begin
      caldat, time, month, day, year, hour, minute, second
      date=long(year)*10000L + long(month)*100L + long(day)
      dateSec=long(hour)*60L*60L + long(minute)*60L + long(second)
    end
    ;Integer: assume this is a month in a monthly climatology
    (timeType eq 2 or timeType eq 3): begin
      date=long(time)
      datesec=replicate(0L, n_elements(time))
    end
  endcase


  ;Write netCDF file
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  fID=ncdf_create(fileName, /clobber)
  
  ; Create required dimensions
  lonDID  = NCDF_DIMDEF(fID, 'lon', long(n_elements(lon)))
  latDID  = NCDF_DIMDEF(fID, 'lat', long(n_elements(lat)))
  timeDID  = NCDF_DIMDEF(fID, 'time', /unlimited)
  
  ; Create a variable to hold the data
  lonID = NCDF_VARDEF(fId,'lon', [LonDID], /float)
  latID = NCDF_VARDEF(fId,'lat', [LatDID], /float)
  dateID = NCDF_VARDEF(fId,'date', [TimeDID], /LONG)
  dateSecID = NCDF_VARDEF(fId,'datesec', [TimeDID], /LONG)
  fluxID = NCDF_VARDEF(fId, 'flux', [LonDID, LatDID, TimeDID], /float)
  
  ; Create some attributes
  NCDF_ATTPUT, fid, /Global, "Title", title
  if keyword_set(author) then begin
    NCDF_ATTPUT, fid, /Global, "Author", author
  endif
  
  NCDF_ATTPUT, fId, lonId, "long_name", "longitude"
  NCDF_ATTPUT, fId, lonId, "units", "degrees_east"
  
  NCDF_ATTPUT, fId, latId, "long_name", "latitude"
  NCDF_ATTPUT, fId, latId, "units", "degrees_north"
  
  NCDF_ATTPUT, fId, fluxId, "long_name", "surface flux"
  NCDF_ATTPUT, fId, fluxId, "units", "kg/m2/s"

  NCDF_ATTPUT, fId, dateID, "long_name", "Date"
  if date[0] gt 10000000L then begin
    NCDF_ATTPUT, fId, dateID, "units", "YYYYMMDD"    
  endif else begin
    NCDF_ATTPUT, fId, dateID, "units", "MM"    
  endelse

  NCDF_ATTPUT, fId, dateSecID, "long_name", "Seconds past 0000 UTC"
  NCDF_ATTPUT, fId, dateSecID, "units", "seconds"
  
  ; Leave definition mode and enter data write mode
  NCDF_CONTROL, fId, /ENDEF
  
  ; Write the data
  NCDF_VARPUT, fId, lonId, lon
  NCDF_VARPUT, fId, latId, lat
  
  NCDF_VARPUT, fId, dateID, date
  NCDF_VARPUT, fId, datesecID, datesec
  NCDF_VARPUT, fId, fluxId, flux
  
  ; Done
  NCDF_CLOSE, fId

  return, 1

end